1. Introduction, Motivation, Relevance, and Objectives

Unemployment remains one of the most visible and impactful indicators of economic well-being, directly shaping the lives of citizens and the priorities of governments. This report provides a comparative analysis of unemployment rates in the United States and Colombia, offering insight into how economic shocks and labor market structures influence short-term employment outcomes.

The United States, with its deep capital markets and relatively flexible labor regulations, often exhibits smooth unemployment cycles responsive to monetary policy. In contrast, Colombia faces more structural unemployment challenges, including higher informality and vulnerability to global commodity price swings.

From a development policy perspective, analyzing and forecasting these trends can inform interventions in education, training, and employment services. This analysis leverages publicly available ILO datasets and advanced time series models to forecast future labor market trajectories and highlight structural contrasts between the two economies.

The main objectives of the report are:

2. Data Description and Preprocessing

The dataset used in this study is sourced from the International Labour Organization (ILO) and contains monthly unemployment data for the United States and Colombia, segmented by sex and age group. To conduct a comparative and consistent analysis, we have constructed a unified time series using the age-group data by summing unemployment counts and calculating a weighted average percentage.

To ensure robust modeling, we identified and corrected anomalies in Colombia’s unemployment rate series. We used the tsclean() function from the forecast package to automatically detect and replace outliers with interpolated values. This method preserves the underlying trend and seasonality while reducing the influence of structural shocks or reporting errors. The cleaned series was then used for all decomposition, model training, and forecasting for Colombia.

Figure: Global Weighted Average vs US and Colombia Unemployment Trends

Figure: Global Weighted Average vs US and Colombia Unemployment Trends

Unemployment rates across the globe have seen significant fluctuations over the decades, often reflecting the impact of global economic crises, policy shifts, and technological change. Historically, an unemployment rate of around 4% is often regarded as full employment, meaning most individuals willing and able to work can find jobs. The global plot shows that many economies experienced sharp spikes during major recessions (e.g., 2008) and more recently during the COVID-19 pandemic in 2020.

However, what’s striking is the rapid recovery of unemployment rates post-COVID in many regions. As reflected in the global and country-level panels, the US shows a pronounced spike in 2020 followed by a fast recovery, thanks to aggressive fiscal and monetary responses. Colombia, while also experiencing a peak, exhibits greater volatility and slower normalization, likely due to structural vulnerabilities such as labor informality and limited social insurance coverage.

These plots collectively set the stage for understanding the differences in labor market resilience between a high-income country and a developing one — an essential motivation for this forecasting exercise.

During initial exploration, Colombia’s unemployment series exhibited irregular seasonal patterns and nonlinear shifts, suggesting the presence of structural outliers. As part of preprocessing, we applied an outlier detection method (IQR-based) and tested both the original and the cleaned series. Visual inspection and model diagnostics showed that the cleaned series led to more stable and interpretable forecasts.

3. Outlier Detection and Pre-Forecast Diagnostics for Colombia

Figure: Original vs Outlier-Removed Series for Colombia (Total Unemployment Percentage)

Figure: Original vs Outlier-Removed Series for Colombia (Total Unemployment Percentage)

The original unemployment series for Colombia (2010–2024) contained several large, abrupt shifts not aligned with typical seasonal or trend patterns. These were likely due to structural shocks—such as policy interventions or labor market disruptions—that introduced high-frequency noise into the series.

Forecast models trained on this unprocessed data displayed:

By applying an IQR-based outlier removal method followed by interpolation, we preserved the core seasonality and trend while stabilizing model input. The cleaned series yielded:

Thus, although both versions were tested, the outlier-removed series was ultimately chosen as the modeling base to ensure robust and interpretable forecasts.

4. Analyzing Decomposition and Stationarity

The decomposed unemployment series for both the US and Colombia (after outlier removal) revealed strong seasonal patterns and structural trends. For the US, the ADF test confirmed non-stationarity (p = 0.5664), while the Mann-Kendall and Seasonal Mann-Kendall tests indicated a significant downward trend and seasonal effects. Similarly, Colombia’s series exhibited non-stationarity (ADF p = 0.63), strong seasonality (Kruskal-Wallis p < 0.001), but no significant seasonal trend (SMK p = 0.136). In both cases, one level of differencing was required to achieve stationarity, and deseasonalized, differenced series were used for robust model training and testing.

The training datasets for both the US and Colombia were constructed by transforming the monthly unemployment percentage into time series objects, excluding the most recent 12 months which were reserved for testing. For the US, data from 2001 to 2023 was used, while Colombia’s training period covered 2010 to 2023 to ensure continuity after handling missing and outlier values. Deseasonalization and differencing were applied where necessary to stabilize the series and meet stationarity assumptions before model fitting.

5 Colombia’s Unemployment Model Fitting:

Replace outliers using tsclean()

## 
##  Augmented Dickey-Fuller Test
## 
## data:  ts_Colombia_total_per_clean
## Dickey-Fuller = -2.6857, Lag order = 6, p-value = 0.2872
## alternative hypothesis: stationary

Step 2: Decompose the cleaned series

Step 3: Fit Time Series Models (using cleaned data only)

## 
##  Ljung-Box test
## 
## data:  Residuals from Seasonal naive method
## Q* = 665.03, df = 24, p-value < 2.2e-16
## 
## Model df: 0.   Total lags used: 24
## Order 1 - 315.6527; Order 84 - 690.3369; Order 168 - 673.258
## Order 1 - 315.6527; Order 42 - 666.0699; Order 84 - 690.3369
## Order 1 - 315.6527; Order 21 - 588.7611; Order 42 - 666.0699
## Order 1 - 315.6527; Order 11 - 488.412; Order 21 - 588.7611
## Order 1 - 315.6527; Order 6 - 409.5131; Order 11 - 488.412
## Order 1 - 315.6527; Order 3 - 353.8308; Order 6 - 409.5131
## Order 1 - 315.6527; Order 2 - 330.5267; Order 3 - 353.8308
## Order 12 - 499.9714

## 
## Model estimated using sma() function: SMA(1)
## Response variable: y
## Distribution used in the estimation: Normal
## Loss function type: MSE; Loss function value: 0.3787
## All coefficients were provided
## Error standard deviation: 0.6172
## Sample size: 168
## Number of estimated parameters: 1
## Number of degrees of freedom: 167
## Information criteria:
##      AIC     AICc      BIC     BICc 
## 315.6286 315.6527 318.7525 318.8143

## 
##  Ljung-Box test
## 
## data:  Residuals
## Q* = 34.924, df = 24, p-value = 0.06954
## 
## Model df: 0.   Total lags used: 24
## 
## Forecast method: Simple exponential smoothing
## 
## Model Information:
## Simple exponential smoothing 
## 
## Call:
## ses(y = deseas_train, h = n_for, holdout = FALSE, silent = FALSE)
## 
##   Smoothing parameters:
##     alpha = 0.8193 
## 
##   Initial states:
##     l = 11.6429 
## 
##   sigma:  0.6107
## 
##      AIC     AICc      BIC 
## 699.1418 699.2882 708.5137 
## 
## Error measures:
##                       ME      RMSE       MAE        MPE     MAPE      MASE
## Training set -0.01120721 0.6071011 0.4487551 -0.2775967 4.258392 0.3883161
##                    ACF1
## Training set 0.01206347
## 
## Forecasts:
##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Jan 2024       10.10021 9.317503 10.88291 8.903164 11.29725
## Feb 2024       10.10021 9.088330 11.11208 8.552674 11.64774
## Mar 2024       10.10021 8.902223 11.29819 8.268048 11.93237
## Apr 2024       10.10021 8.741371 11.45904 8.022046 12.17837
## May 2024       10.10021 8.597641 11.60277 7.802230 12.39818
## Jun 2024       10.10021 8.466507 11.73391 7.601678 12.59874
## Jul 2024       10.10021 8.345144 11.85527 7.416070 12.78434
## Aug 2024       10.10021 8.231647 11.96877 7.242491 12.95792
## Sep 2024       10.10021 8.124660 12.07575 7.078868 13.12155
## Oct 2024       10.10021 8.023176 12.17724 6.923663 13.27675
## Nov 2024       10.10021 7.926425 12.27399 6.775695 13.42472
## Dec 2024       10.10021 7.833801 12.36661 6.634038 13.56638

## 
##  Ljung-Box test
## 
## data:  Residuals from Simple exponential smoothing
## Q* = 33.538, df = 24, p-value = 0.09322
## 
## Model df: 0.   Total lags used: 24
## Series: ts_train 
## ARIMA(2,0,0)(0,1,1)[12] 
## 
## Coefficients:
##          ar1     ar2     sma1
##       0.8197  0.1305  -0.8273
## s.e.  0.0801  0.0803   0.0789
## 
## sigma^2 = 0.4384:  log likelihood = -162.86
## AIC=333.72   AICc=333.99   BIC=345.92

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,0,0)(0,1,1)[12]
## Q* = 37.599, df = 21, p-value = 0.01436
## 
## Model df: 3.   Total lags used: 24
## Series: deseas_train 
## ARIMA(1,1,0) 
## 
## Coefficients:
##           ar1
##       -0.1525
## s.e.   0.0771
## 
## sigma^2 = 0.3744:  log likelihood = -154.45
## AIC=312.89   AICc=312.96   BIC=319.13

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,1,0)
## Q* = 34.722, df = 23, p-value = 0.05541
## 
## Model df: 1.   Total lags used: 24

## 
##  Ljung-Box test
## 
## data:  Residuals from STL +  ETS(M,N,N)
## Q* = 53.393, df = 24, p-value = 0.0005132
## 
## Model df: 0.   Total lags used: 24

## 
##  Ljung-Box test
## 
## data:  Residuals from Regression with ARIMA(2,1,0) errors
## Q* = 126.3, df = 22, p-value < 2.2e-16
## 
## Model df: 2.   Total lags used: 24

## 
##  Ljung-Box test
## 
## data:  Residuals from TBATS
## Q* = 27.782, df = 24, p-value = 0.2694
## 
## Model df: 0.   Total lags used: 24

## 
##  Ljung-Box test
## 
## data:  Residuals from NNAR(3,5)
## Q* = 39.214, df = 24, p-value = 0.02593
## 
## Model df: 0.   Total lags used: 24

## 
##  Ljung-Box test
## 
## data:  Residuals
## Q* = 38.377, df = 24, p-value = 0.03171
## 
## Model df: 0.   Total lags used: 24

## 
##  Ljung-Box test
## 
## data:  Residuals from StructTS
## Q* = 186.02, df = 24, p-value < 2.2e-16
## 
## Model df: 0.   Total lags used: 24

Step 4: Accuracy Comparison for All Models (Cleaned Data)

## ✅ Best model based on Average(RMSE, MAPE): NNETAR
Forecast Accuracy on Cleaned Data (Average Based on RMSE and MAPE Only)
ME RMSE MAE MPE MAPE ACF1 Theil’s U Average
SNAIVE 0.042 0.650 0.525 0.214 5.394 0.210 1.080 3.022
SMA -0.583 1.300 1.108 -7.582 12.159 0.706 2.501 6.729
SES -0.467 1.252 1.050 -6.362 11.454 0.706 2.352 6.353
SARIMA -0.635 0.857 0.711 -7.239 7.954 0.574 1.708 4.405
ARIMA -0.500 1.267 1.068 -6.712 11.668 0.705 2.394 6.467
ETS -0.512 0.811 0.661 -5.978 7.378 0.613 1.609 4.095
ARIMA + Fourier 0.014 0.566 0.464 -0.439 4.871 0.479 1.002 2.718
TBATS -0.229 0.510 0.435 -2.761 4.742 0.392 0.977 2.626
NNETAR 0.083 0.417 0.312 0.685 3.279 -0.052 0.747 1.848
SSES -0.086 0.462 0.387 -1.288 4.175 0.349 0.872 2.319
BSM 0.689 0.868 0.718 6.783 7.108 0.562 1.446 3.988
Average of 3 -0.077 0.432 0.352 -1.121 3.788 0.226 0.806 2.110

Step 5: Forecast for 2025 using NNETAR, SSES, and TBATS

## [1] 9.03823